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We analyze the Gaussian approximation as a method to obtain the first and second moments of a stochastic 
process described by a master equation. We justify the use of this approximation with ideas coming from van 
Kampen's expansion approach (the fact that the probability distribution is Gaussian at first order). We analyze 
the scaling of the error with a large parameter of the system and compare it with van Kampen's method. Our 
theoretical analysis and the study of several examples shows that the Gaussian approximation turns out to be 
more accurate. This could be specially important for problems involving stochastic processes in systems with a 
small number of particles. 

PACS numbers: 



INTRODUCTION 

Master equations are a convenient tool to treat stochastic Markov processes!]! 01 . In some cases, they offer an alternative 
approach to the Chapman-Kolmogorov equation and have been used extensively in discrete-jumps, or birth-death, processes, 
such as chemical reactions (including those happening inside a cell), population dynamics or other ecology problems[3], opinion 
formation and cultural transmission in the field of sociophysics[4], etc. In all these cases, it is important to consider that the 
population number (whether molecules, individuals, agents, etc.) might not be very large (maybe ranging in the tens or hundreds) 
and the fluctuations, that typically scale as the square root of the inverse of this number, can not be considered as negligible. It is 
therefore, of the greatest importance to derive evolution equations for the average behavior and the fluctuations. The important 
work by van Kampen[ l] offers a systematic way of deriving these equations from an expansion of the master equation in a 
parameter i2, typically the system volume. The £2-expansion is mostly used in its lowest order form, in which one can prove 
that the error in the average value, the second moment and the fluctuations (the variance), scale at most as OP, Q} and H 1 / 2 , 
respectively. The van Kampen il-expansion, furthermore, shows that, at this lowest order, the fluctuations follow a Gaussian 
distribution. In this paper, we take this result of van Kampen's theory and, considering from the very beginning that fluctuations 
are Gaussian, we derive a closed system of equations for the average value and the second moment. This Gaussian closure of the 
hierarchy of moments turns out to be more accurate that the il-expansion as the above-mentioned errors scale at most as QT l l 2 , 
£l l l 2 and H 1 / 2 , respectively. Furthermore, the Gaussian closure scheme is very simple to carry on in practice and can be easily 
generalized to systems described by more than one variable. 

The paper is organized as follows: In the following section, we will briefly review the £2-expansion and derive the main 
equations for the Gaussian closure approximation. The errors of both methods are discussed in section „ In sections „ and „ we 
will give examples of the application of the method in the cases of a binary chemical reaction and an autocatalytic reaction. 
The results of these two examples confirm the error-analysis performed before. For both processes we compare with the results 
coming from the exact solution of the master equation in the stationary regime (derived in the appendix for the binary chemical 
reaction), and the results of numerical simulations using the Gillespie algorithm in the time-dependent evolution. In section n we 
present an application to a recently introduced model for opinion formation which requires two variables for its full description. 
Finally, in section„we end with a brief summary of the work. 



FORMULATION 



Let P(n,t) be the probability that at time t the population number takes the value n. We consider that it evolves according to 
a general master equation of the form: 

= l)[C k (n;£l)P(n,t)], (1) 

where E is the linear step operator such that E k [fin)] = f(n + k) and k runs over the integer numbers. Besides n, the coefficients 
Q (n;i2) depend on Q., which is a large parameter of the system (typically the system volume). We consider that these functions 
are polynomials or can be expanded in power series of n as Q(«;i2) = £ fl C^(n)n" where the coefficients Ci{£l) scale as 

C|(fl) = ( c a k(j + c a k j£2~' +c£ 2 £2~ 2 + ...). Master equations of this form appear in the description of chemical reactions 
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JH, ecological systems ^ and opinion dynamics 0, among many other cases. More specific examples will be considered in 
the next sections. 

In a seminal work, van Kampen[ 1] has given a way of finding an approximate solution of Eq. (Q]). The approximation is based 
upon the splitting of the variable n using the ansatz n = £20 (r) +£22^, where 0(f) ~ <9(£2°) is a function of time accounting 
for the deterministic part of n and t, ~ (9(£2°) corresponds to the fluctuations. Changing variables from n to | in Eq. {T), and 
expanding in powers of £2 one obtains a Fokker-Planck equation for the probability distribution II(§ ,t) of the new variable % : 
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dt 
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where the macroscopic variable satisfies 
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From Eq.© we obtain the first and second moments of the fluctuations: 
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As proven by van Kampen, the solution of the Fokker-Planck equation (|2j is a Gaussian distribution. Therefore, the £2- 
expansion method tells us that, up to corrections of order £2 2, the fluctuations of the variable n follow a Gaussian distribution. 
It suffices, then, to know the first and second moments of this distribution. Our intention is to use from the very beginning the 
Gaussian property in order to obtain a closed system of equations for the first two moments (n) and (n 2 ). 

From ((TJ we get the (exact) equations for these first two moments, as: 

tH ""' :£<*(*-B)Q(n;n)). (6) 



d{n) 



dt 



= -£<£Q(n;£2)>, 

k 



dt 



After substitution of the series expansion Q(«;£2) = £„Q'(£2)n a in the right hand side of these equations, one obtains higher 
order moments («'") for m > 3. The Gaussian closure replaces these higher order moments with the expressions (n m )c that hold 
in the case of a Gaussian distribution, i.e. (n)c = (»), (" 2 )g = (« 2 ) an d 



(n m )c 



k=l 



(2k- 1 



\\\(n\ m ~ 2k 



[(n 2 



(7) 



for m>3. The first moments are explicitly shown in table U 



Moment 


Gaussian approximation 


<« 3 > 


3<n 2 )(rc)-2(n) 3 


{n 4 } 


3{n 2 ) 2 -2{n) 4 


{n 5 } 


\5{n 2 ) 2 {n) -20(n 2 )(n} 3 +6(n} 5 


(nO) 


15(rc 2 ) 3 -30<rc 2 )<rc) 4 + 45(rc) 6 


(n 2 n 2 ) 


[n\)(n 2 ) + 2{n l ){n l n 2 ) -2(n x ) 2 (n 2 ) 


(nfnf) 


(n 2 )(n 2 }+2(n l n 2 } 2 -2(nl} 2 (n2) 2 


(»1«2> 


3{n 2 ){n x n 2 ) -2(wl) 3 (n2) 


(n\n 2 ) 


6{n l n 2 ) 2 {n x )+6{n\)\n2) 2 + 6{n i n2){n 2 ){{n x ) 2 -2{n l ) 2 
-6{n 2 ){n 2 ) 2 {n l )+3{n 2 ){n 2 )(n l )-2{n l ) 3 {n 2 ) 



TABLE I: Gaussian moments 



The van Kampen ansatz n = £20 (f ) + £22 1 allows us to find the error of this approximation. It follows that: 



£2'"- 1 



E 



£2 1 -'/ 2 0"'- / (4') 
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In the Gaussian approximation, the first three terms of the sum, / =0,1,2 are exact and the term / = 3 scales as Q. 2 , or: 



^^+o(a-^). ( 9) 



ft'"- 1 £1 

If we use this result in each of the terms of Eq.© and Q.(h;£2) = Q l ~ a {c k ' + 0(£2 -1 )) we obtain 

d(n 



dt 

with g\ = — \\ {kCk{n;Q)) G . Similarly, one finds 



((»>,(» 2 )) + 0(ft- 1 / 2 ), (10) 



k 

din 2 ) 



dt 



:g 2 ((»>,(« 2 ))+0(ft 1 /2 )j (U) 



with g 2 = Y^{Kk-n)C k {n-a)) G . 

k 

This Gaussian approximation scheme (or equivalently, finding a hierarchy of equations for the cumulants and neglecting those 
of order greater than two) has been used many times in the literature in different contexts 0,01 • We will show in the next section 
that the direct use of Eqs. d 1 01 1 lb has a smaller error that the use of Eqs. (OH). Before showing this, we will generalize this 
procedure for the case of two-variable problems. Let us consider a master equation of the following form: 

dP( ni n 2 ,t) = ^ {E h E h_ l) [c kuk2 {n u n 2 ;a)P{n h n 2 ,t)] . (12) 

The evolution equations for the first, second order moments and the correlations are: 

dlrii) 



dt 

d(nf) 
dt 

d{ri\n 2 ) 
dt 



= - (*Ai,fe(»i)»2;^))) (13) 
k u k 2 

= L {ki{ki-ni)C klM {nun 2 ;Q)), (14) 
ki ,k 2 

= {{k\k 2 -k 2 n\ -k\n 2 )C kukl (n\,n 2 ;Cl)) , (15) 



(; = 1,2). Again, the Gaussian closure consists in replacing (n'"'n™ 2 ) by the expression (h^'h^g tnat holds assuming that the 
joint distribution P(n\,n 2l t) is Gaussian. This can be computed using Wick's theoremfl. In table (HJ) we write the expression of 
some of the terms. 



ERROR OF THE METHOD 



We now calculate the error of the Gaussian approximation and compare it with the one of the H-expansion. In Eqs. ( fTOlfTTb 
we have shown that the errors we introduce in the equations for the moments when performing the Gaussian approximation are 
of order 0{£l~ x l 2 ) for (n) and 0{£l l l 2 ) for (n 2 ). The Gaussian approximation scheme proceeds by considering approximations 
Hi(t), Hi{t) to the true moments (n(t)), (n 2 (t)). These approximations are defined as the solution of the evolution equations 
(fToTTTT ): 

^-=glGlll,M2), ^=#2(Ml.M2)- (16) 

Defining the errors E\ , E 2 as: (n) = fJ-i + £\ , (n 2 ) = [i 2 + e 2 ; expanding in first order in E\ and e 2 , and using equations ( fToTTTTb 
and (fTol l we get: 

de, = d gl (^ 2 ) d gl fa,K) £2 + Q(n -i/ 2)) (17) 

dt djli djj, 2 

de, = dg^al + £2 + 1/2 

dt dfl\ djl 2 
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Taking into account that jii ,g\ ~ (9(£2), 1^2, g2 ~ <9(£2 2 ), we have: 

^ = 0(i2°)e 1 + 0(i2- 1 )e 2 + 0(i2- 1/2 ), (19) 

^ = 0{a)e x +0{QP)e 2 + 0{D}l 2 ). (20) 

If we set £i ~ <9(£2"), £2 ~ 0(Q. h ), and the initial conditions are known, so that initially £1 = £2 = 0, equations dT9b , ( 120b imply 
that a < —1/2 and b < 1 /2, a scaling respected during the time evolution. 

In conclusion, solving equations ( fTOlfTTl l. we get (n) and (n 2 ) with errors of order £1 = (9(£2~ 1//2 ) and £2 = 0{£l x l 2 ), or 
smaller. Using the equations ((3]|5]l of first order van Kampen's expansion the error is of higher order in both cases: 0(£2°) for ( n) 
and 0{Q}) for (n 2 ). However, for the variance, a 2 = (n 2 ) — (n) 2 , both approximations have an error of order <9(n'/ 2 ). We will 
show in the next sections that the Gaussian approximation has the extra advantage that it is easier to derive for many problems 
of practical interest. 

One might be tempted to go to higher order schemes, where one neglects all the cumulants of order greater than m with m > 2, 
and in this way obtain a closed set of equations for the first m moments. For example, if we neglect all the cumulants of order 
greater than 3, applying the same analysis as before, it is possible to derive that the errors in the first, second and third moments 
are of order (9 (£2 ~ 1 , £2° , £2 1 ) , respectively. 

A word of caution is needed here. When truncating beyond the second cumulant, it is not ensured that the resulting probability 
distribution is positive definite ifioll . This means that one could get from such an scheme inconsistent results, e.g. a negative 
variance. Nevertheless, according to our analysis, the importance of these spurious results would decrease with £2 as indicated, 
so one can still get useful results from higher order schemes. 

In the following sections we will compare the Gaussian approximation presented here with the first order £2-expansion in 
some specific examples. 



BINARY REACTION A + B-±C 

Chemical reactions are suitable processes for a stochastic description. The stochastic approach is specially necessary when 
the number of molecules considered is small, as it is the frequently addressed case of chemical reactions inside a cell, because 
in this situation fluctuations can be very important. 

K 

We consider the general process A + B ^± C, limited by reaction. This means that any two particles A and B have the same 

probability of reaction. Denoting by A(t) and B(t), respectively, the number of molecules of the A and B substances, the rate for 
the A+B — > C reaction is -]^A(t)B(t). For the reverse reaction, it is assumed that C has a constant concentration, and hence the 
rate is ©£2. In these expressions £2 is proportional to the total volume accessible. Since B(t) — A(t) = A is a constant, one only 
needs to consider one variable, for example, the number of A molecules at time t . Let us denote by P(n,t) the probability that 
there are n A-molecules at time t . The master equation describing the process is: 

dP(n t) K 

dt = ^[( n + 1 )( A + n + 1 ) P ( n + 1 '0-»(» + A)P(n,f)] + ^[^(n-l,0-^(n 5 0] I (21) 

which is the basis of the subsequent analysis. Note that this equation can be written in the form (fTJ setting C\ (n;£l) = ^n(n + 
A),C_i(n;£l) = fl)Q. 

In the irreversible case, CO = 0, this master equation can be solved exactly using the generating function technique. In the 
general case, CO ^ 0, an exact solution can also be found for the stationary state dP g"/^ = 0. Details of the calculation are given 
in the appendix. We will compare the results obtained from the Gaussian approximation and the first order £2-expansion with 
the exact results, when available. 

The equations for the first two moments, using (O, are: 

d(n) 



dt £2 

din 2 



= -p:((» 2 )+A(n»+£2(0, (22) 



= -(-2(n J ) + (l-2A)(« 2 )+A(«))-2£2o)(n)+£2co. (23) 



dt £2 

Using the Gaussian approximation, the evolution equations for the moments are: 



= "^(M2+AMi)+£20), (24) 
A = ^(4Ati 3 "6^2Ati + (l-2A)Ai 2 + AAii)+2£2a)Aii+£2c(J. (25) 
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And the first order £2-expansion gives: 



d(j) 
dt 

d(n) 

dt 

d(n 2 ) 
dt 



-K(j)(<p + 8) + <a, 
KD.(j) 2 -K(d + 2(j))(n) +Q.CO, 

~2K(2<j) + 8)(n 2 ) +Q[jc0(0 + 5)(l-2(n)) 
2(fc£20 2 (2</> + 8) + ©((«) + 1 + cotj)))] , 



(26) 
(27) 

(28) 



where 8 = ^. 



We compare the two approximations in the time-dependent case with results obtained by averaging over single realizations of 
the process, obtained numerically using the Gillespie algorithm! 3]. In the next figures we compare the exact results with those 
obtained from the Gaussian approximation (computed by numerical integration of equations l24l l25l l and £2-expansion (equations 
I26T28I 




— Gaussian 
^ ■ £2- Expansion 

Gillespie 
Exact stationary 




6 9 
t 

FIG. 1: (n(t)}, (n 2 (t)} and o 2 (t) for the binary reaction A + B ^± C with parameters K= 1, to = 1, £2 = 10 and initial conditions n(0) = 100, 

5 = 1. For the first two moments the Gaussian approximation (solid) is very close to the results obtained with the Gillespie algorithm (dot- 
dashed, obtained averaging over one million realizations) and the exact stationary value (thin line), while I st order fl-expansion (dashed) gives 
clearly different values. For a 2 , the £2-expansion gives more accurate results but both approximations differ from the exact values. 



2 2 

Error in <n> Error in <n > Error in a 




1 10 100 1000 10000 1 10 100 1000 10000 1 10 100 1000 10000 
Q. Q. Q. 



FIG. 2: Error in (n) , (n 2 ) and a 2 in the stationary state in the same case than in FigQ] The straight thin lines are fits to the data and have slope 
— 1, or 1. For the Gaussian approximation (solid), the errors in ((n), {n 2 ), rj 2 ) scale as (n _1 , £2°, ft ). For the Q. -expansion (dashed), the 
errors scale as (£2°, £l l , £2°). 

Figure (Q]i shows that the Gaussian approximation reproduces better the exact results for the first two moments; for the 
variance, the £2-expansion gives more accurate results but both approximations differ from the exact values. Figure (|2|i shows 
that the errors in the stationary state, coming from the Gaussian approximation for the mean value, the second moment and the 
variance scale as (£2 _1 , £2°, £2°), respectively, while the errors of the £2-expansion at first order scale as (£2°, £2', £2°). This 
scaling is consistent with the previous analysis, as the exponents of the errors are smaller than the obtained bounds. 



AUTOCATALYTIC REACTION A — X, 2X -» B 



The master equation describing this process is[Q]]: 
dP(n,t) 



k J 



dt 



Q.fak[P(n-l,t)-P(n,t)]+-[(n + 2)(n + l)P(n + 2,t)-n(n-l)P(n,t)], 



(29) 



where the concentration of A particles is consider to be constant with a value fa. This equation if of the form (Q]) withC_i(n;£2) = 
£lkfa,C2(n;Q.) + j>n(n — 1). The general solution for this equation is not known, but the stationary solution P st (n) can be 
obtained using the generating function technique! 1]. The exact equations for the first moments are: 



_ = £l kfa+2k — -2k—, 
^ = ii^(2(„) + l)-^(4(„ 3 ) 



(n 2 )+4(n)). 



Performing the Gaussian approximation, we get: 



= Clkfa+lk'^-lk'^, 
dt m £2 £2 

- Q^ a (2^ 1 + 1)--(12^ 2 Mi-8Mi-8M2 + 4mi)- 



While first order £2-expansion approach leads to: 

— = MA -2k 1 it) 2 , 
at 



din) 
dt 
d(n 2 ) 
dt 



= £2(/c0 A +2/c> 2 )-4/c>(n), 

= -8/c>(« 2 > +Q.(2kfa+4-k'<t> 2 )(n) +Q.(kfa +4/c> 2 ) 



(30) 
(31) 

(32) 
(33) 

(34) 
(35) 
(36) 



In the next figures we show the results obtained with the Gaussian approximation (computed by numerical integration of 
equations [32ll33l , £2-expansion (equations [3411361 ). the Gillespie algorithm, and the exact stationary solution. 



2 

<n > 



10.1 




1 1 1 1 


f 










1 



10 




FIG. 3: (n(t)), (n(t) 2 ) and cx 2 (f) for the autocatalytic reaction A -» X, 2X -> B withkfa = 1, Id = 1/2, D. = 10 and initial condition n(0) = 0. 
For the first two moments the Gaussian approximation (solid) is very close to the results coming from the Gillespie algorithm (dot-dashed) and 
the exact value in the stationary case (thin line) whereas the fl-expansion result (dashed) is clearly different, although for o" 2 the f2-expansion 
provides more accurate results. 



As in the previous example, we see that the Gaussian approximation fits better the evolution of the moments, but the variance 
is somehow better approximated by the first order £2-expansion. In figure (|4|i we show the errors in the stationary state for the 
two approximations as a function of £2. We see that the errors in ((h), (« 2 ), a 2 ) decay as (£2 , £2 ,£2°) for the Gaussian 
approximation, while the first-order £2-expansion leads to errors that scale as (£2°, £2 1 ,£2°). Again, this scaling is consistent with 
the analysis of the approximations performed. 



7 



2 2 
Error in <n> Error in <n > Error in a 
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FIG. 4: Error in (n),(n 2 ) and o" 2 in the stationary state as a function of £2 in the same case than in Fig(3] The thin lines have slope — 1, or 
1. For the Gaussian approximation (solid), the errors in ((n), (n 2 ), a 2 ) scale (asymptotically) as , i2~', £1°). For the £2-expansion, the 
errors scale as fi 1 , OP). 



OPINION FORMATION 



In the last few years there has been a growing interest in the application of methods and techniques coming from statistical 
physics to the study of complex phenomena in fields traditionally far from physics research, particularly in biology, medicine, 
information technology or social systems. In particular the application of the physical approach to social phenomena has been 
discussed in several reviews fi , 11 , 12 ] . As an example of the use of master equations in this field, we mention a recent paper 
in which the process of opinion formation in a society is modeled as follows: Society is divided in two parties, A and B, plus 
an "intermediate" group of undecided agents I. The supporters of A and B do not interact among them, but only through their 
interaction with the group I, convincing one of its members with a given probability. In addition there is a nonzero probability 
of a spontaneous change of opinion from I to the other two parties and vice-versa. More specifically, if n A ^ is the number of 
supporters of party A(B), «/ is the number of undecided agents and Q. is the total number of individuals, the possible transitions 
are: 

• /, occurring with a rate a\n A , 

A, occurring with a rate a2«/, 

• /, occurring with a rate a^n B , 

B, occurring with a rate o^n/, 
2A, occurring with a rate ^n A rij, 
2B, occurring with a rate ^n B ni. 

«a + n B + «/) is fixed, there are only two independent variables, say n A and n B . The 



spontaneous change A 
spontaneous change / 
spontaneous change B 
spontaneous change / 
convincing rule A + 1 

convincing rule B + I - 
As the total number of individuals (Q. 
master equation of the process is: 



—P{n A ,n B ,t) = ai(n A + l)P(n A + l,n B ,t) + a3(nB + l)P(n A ,n B + l,t) 
+a2{£l-n A -n B +l)P(n A - I ,n B ,t) + a^Sl - n A - n B + l)P{n A ,n B - l,f 
+{Q.-n A -n B + \) 



(37) 



^■(n A -l)P(n A - l,n B ,t) + ^(n B - l)P (n A ,n B - l,t) 



a\n A + 0}Hb + (0C2 + «4 J (« - «a - n B ) H — (Ll-n A -n B ) 



P(nA,n B ,t) 



We note that this master equation can be written in the general form (fT2] i by setting C\ = <X\n A , Co,i = d^n B , C_i.o = 

(Q.-n A -n B )(a 2 + andCo,-i = {Q.-n A -n B ){04 + ffns). 

An exact solution of this master equation is not known. In the following, we will apply to this problem the Gaussian approxi- 
mation scheme and compare it with the results of the H-expansion. The exact equations for the first moments are: 
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d(n A (t)) 
dt 

d(n B (t)) 



dt 

nil 



d{n\{t)) 



dt 



d{n 2 B {t)) 
dt 



d{n A (t)n B {t)) 
dt 



- (at + 052 - ft) (n A ) + «2 (Q - (n B ) ) - ^ - § («a«b) , 

- (0£ 3 + «4 - ft)(«s) + 054 (£2 - (%) ) - ^ - ^ («A«b) . 

(a! + a 2 (2£2 - 1) + ft)(« A ) + a 2 (£2 - (n B ))- 2{a x + a 2 - ft + ^)<«1) 
-(2a 2 + ^-)(«a«b) - - ^-M»b), 

(a 3 + a 4 (2£2 - 1) + ft)(/i B ) + a 4 (£2 - ("a» - 2(a 3 + a 4 - ft + 

^ , ft s , , 2ft , 3 2ft 2 \ 
- (2a 4 + — ) {n A n B ) --Q\ n B)--Q ( n A n Bl > 

-(ai + a 2 + a 3 + a 4 - ft - p 2 )(n A n B ) + a 2 {D.(n B ) - (n|)) 



+a 4 (£2(« A )-(« A ))- 



ft+ft 
£2 



{{n 2 A n B ) + {n A n 2 B )). 



(38) 
(39) 

(40) 

(41) 

(42) 



Denoting by Ai,A 2 ,Z?i,Z? 2 ,C the Gaussian approximations to the moments (n A ), (« A ), («b), («|) and the correlation (n A n B ) 
respectively, and using the results in table H we obtain: 



dAi 

dt 
dB { 

dt 
dA 2 
dt 



dB 2 
dt 



dC 
dt 



-(oci + a 2 - ft)Aj + a 2 £2- a 2 B! - ^-A 2 - ^-C, 

ft ft 

-(a 3 + a 4 - ft)Bi + a 4 £2 - avh - gfi 2 - gc, 

(oc 1 + oc 2 (2£2-l)+ft)A 1 + oc2(£2-Bi)-2(ai + of2-ft + ^)A 2 

- (2a 2 + ^ )C - ^ (3A iA 2 - 2A^) - ^ (A 2 B i + 2AiC - 2A?2?i) , 

ft 

( a 3 + a 4 (2£2 - 1 ) + ft )B , + a 4 (£2 - A , ) - 2 (a 3 + a 4 - ft + g )B 2 
-(2a 4 + ^)C- ^(3Bi2?2-2B?)- ^(BaAi+2»iC-2flfAi), 



-(ai + a 2 + a 3 + a 4 - ft - ft)C + a 2 (£2fi! - B 2 ) 
+a 4 (£2A! -A 2 ) 



(43) 
(44) 

(45) 

(46) 
(47) 



^ " ! * [BiA 2 +B 2 A 1 +2(A 1 +A 2 )C-2A?A 2 -2B?B 2 ] 



£2 

In van Kampen's expansion method, we define 0a(b)i4a(b) sucn tnat n A(B) — ^0a(b) + ^ 1 ^ 2< =a(b)- 
The equations for the macroscopic components are jg]: 

= -ai0 A + [O2 + ft0A](l-0A-0B), 

dt 

= -a 3 <fe + [a 4 + ft<fe](l-</> A -<fe), 



(48) 
(49) 
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and for the fluctuations: 



dt 



dt 

£2 



dt 



= -[a 1 + a 2 + /3i(20 A + 0iO-M4A)-(a2 + /3i0A)(&>, (50) 
= -[o& + 054 + ft(2^B + ^ A )-ft]<&)-(a» + ft*B)<6i). (51) 

= -2ai<#> -2((fe+M0«^) + (^6,)) + 2J5i<^)(l - fc- fe) 

+ai& + (Gk + j8i&)(l-&-fo), (52) 

- -2a 3 (ll) -2(a4 + j8 2 B )((§l) + (^))+2j3 2 ^ 2 )(l -fc- <fe) 

+a 3 <fe + (o^ + j3 2 <fe)(l-0 A -<fe), (53) 

= -(ai + osX^) - (c& + /3ito)«&§») + (<^ 2 )) - (a 4 + /3 2 ^)«^^> + 

+(1-0a-^)(J3i+j8 2 )(^ b ). (54) 



From those we can recover the original variables n A /g\(t), 

In the next figures we compare the results coming from both approximations (obtained by numerical integration of the previous 
equations) and from simulations of the process using the Gillespie algorithm, for some representative values of the parameters 
and initial conditions. Again, the Gaussian approximation reproduces better the values for the average and the second moment 
whereas in this case both methods perform very similarly for the fluctuations and correlation. 




FIG. 5: (« A (f)}> ( w l(0)> °A( f ) an d °ab(0 — { n A n B} — ( n A)( n B) f° r the opinion formation model of reference @], for a; = /3; = 1,£2 = 10, 
and initial conditions n^(0) = 0, «b(0) = £2. For the average (n A (f)}, the Gaussian approximation (solid) follows very accurately the Gillespie 
simulation results (dot-dashed), whereas the fl-expansion (dashed) differs clearly. For the second moment (ha (t) 2 ) the Gaussian approximation 
performs clearly better as well, while for the variance cr|(f) and correlations cr| B (?), the Gaussian approximation and the ^-expansion give 
very similar results, although both are far from the simulation data. 



CONCLUSIONS 



In this paper, we have given explicit expressions for the equations for the first and second moments of a stochastic process 
defined by a general class of master equations using the Gaussian approximation closure. The approach is motivated by van 
Kampen's £2-expansion result that, at lowest order, the fluctuations are Gaussian. We have shown that the Gaussian closure is 
simple to perform and leads to errors in the average value, the second moment and the fluctuations (the variance), that scale at 
most as (£1~'/ 2 , H 1 / 2 , £2 1//2 ), respectively. This is to be compared with the £2-expansion result in which the respective errors scale 
at most as (OP, Q} , H 1 / 2 ). Therefore, the Gaussian approximation is more accurate, which turns out to be important, specially 
for small values of £2. We have checked these results by comparing the performance of the two methods in three examples: (i) a 
binary chemical reaction, (ii) an autocatalytic reaction and (iii) a model for opinion formation. In all cases studied, the Gaussian 
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closure has given a better approximation to the average and the second moment, although the i2-expansion, due to a cancellation 
of errors, yields a somehow smaller numerical error in the variance. In general, the Gaussian closure scheme is very simple to 
carry on in practice and this simplicity and the improvement of the predictive power is more apparent in many-variable systems. 
We believe that this method can be usefully applied to the study of other problems of recent interest in the literature involving 
stochastic processes in systems with a small number of particles. 



APPENDIX: REACTION-LIMITED PROCESS 

We now find the solution of the master equation ( f2TT > in the equilibrium state for the general case, and the full dynamical 
solution for the irreversible case (0 = 0. Without loss of generality, let us rescale t — >• Kt/Cl and (0 —> (OQ. 2 / K to get the simpler 
equation: 



dP(n,t) 
dt 



(n + l)(A+n + l)P(n + l,t)-n(n + A)P(n,t) + 0)[P(ri-l,t)-P(ri > t)]. 



(55) 



Furthermore, only the case A > needs to be considered. If A < the change nf = n — A leaves invariant the previous equation 
provided that we make the identification P(n,t) —> P(n + A,t). This means that the solutions in both cases are related by 
P(n,t;A) =P(n-A,t;-A). 
The generating function 



f(s,t)=Y,P(n,ty, 

n=0 



(56) 



satisfies the partial differential equation: 



ff _4 + (l+A)/-tD/ 
as 1 - as 



(57) 



Let us first discuss the equilibrium solution in the general case. 



The equilibrium solution 

By setting ~£ = one gets the differential equation: 

d 2 f , sdf 

s -4 + (l+A )-J—(Df=Q. (58) 
as as 

The solution around the singular regular point 5 = can be found by the Frobenius method as a power series Lr=o fl « s " +V - The 
regular solution satisfying the boundary condition f(s = 1) = 1 is ll 1 3fl : 

f(s)= . J , (59) 

I A (2^ (OS) 

and the equilibrium probabilities, rescaling back to the original parameters, are: 

ctf^r* , (60) 

I h [2Q.^ajKjn\{n+A)\ 

from where the first two moments can be computed as: 

7 A+ i (2Q.^(oJk) I h+l (2Q.^(ojK] 

• '-Q-y/cyic, (n 2 )=Q. 2 (o/K ^ ^-AD.^/a/K. (61) 
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The time-dependent solution 



We now study how the system relaxes towards equilibrium. We will restrict ourselves to the irreversible case CO = 0. This 
corresponds to the process A + B — > 0, inert. The partial differential equation d57l i can be solved by the technique of separation 
of variables by trying solutions of the form f(s,t) — f\ (s) f2(t). This leads to the pair of ordinary differential equations: 

s(l-s)f[' + (l-s)(l+A)f[ + X 2 f 1 - 0, (62) 
/2 + AV2 = 0, (63) 

being A the constant arising from the method of separation of variables. The solution of the time dependent function is e ' 
and the solution of the ^-function is the hypergeometric function II 1411 F(— /ii,/i2',A+ l;s). The explicit series is: 

(a)„ is the Pochhammer's symbol: (a) n — ^fe^p, or (a)o = 1, (a) n = a{a + 1) . . . (a + n — 1) for n > 0, and we have introduced 



-A + VA 2 +4A 2 A+VA 2 +4A 2 
Mi = 2 ' M2 = 2 ' (65) 

The solution for the function f(s,t) is obtained by linear combination of the elementary solutions found above: 

f(s,t) =Y J CxF(-^ l ^ 2 ;A+Us)c- x2 ' . (66) 
X 



This function is, in general, an infinite series on the variable s. In fact the coefficients, according to ( 1561 ) are nothing but the 
time-dependent probabilities. However, in this irreversible case, the probability of having more A-molecules that the initial 
number at t = 0, say M, has to be zero. Therefore the series must be truncated after the power s M . This implies that in the 
previous expression only hypergeometric functions that represent a polynomial in s can be accepted. This is achieved by forcing 
Hi = k = 0, 1,2 . . . ,M, since the series d64b becomes then a polynomial of degree k. The condition jii — k is equivalent to the 
parameter A adopting one of the possible values = y/k(k + A). Finally, noticing that /I2 — Mi = A, the solution can be written 
as: 

M k 

/M = L L C k (A,M) e - k ^'B n , k (A)s n . (67) 

k=0n=0 

The notation emphasizes that Q depends both on A and M but B„ j depends only on A: 

s , ii(a) = HM^. (68) 

n!(A+l)„ 

All that remains is to impose the initial condition. We start with M A-molecules at time t = 0, such that f(s.t — 0) = s M . This 
implies that the coefficients Q must satisfy: 

M 

52 = <S„,m, (69) 

for n = 0, 1, . . . ,M. The solution starts by finding first Cm = 1 /Bm,m and then proceeds backwards to find Cm-i,Cm-2, ■ • ■ ,Co in 
a recursive manner. After some lengthy algebra, the result is: 

ky ' ' V ; £ + A A! (M + A+iy V 7 

(in the case A = k = the correct interpretation of the undetermined expression is Co = 1). Going back to the original time 
variable, we now give the expression for the probabilities: 

M 

P{n,t) = £ C k (A,M)B n .k(&)e- k{k+A)Kt/a - (71) 

k=n 
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k 

The normalization condition Ym=o^n(t) = 1 is verified with the help of the relation £n=o^n.* = (^o- The relation ^ nB,,^ = 
, A! 

(—l) k—— (the indetermination arising when A = k = must be resolved as 0) helps to find the average of the number of 
particles: 

<n(0> = £ (2fc + A) i ^lll 1 fc - kik+A) ' a/R (72) 
The second moment (n(t) 2 ) can be found with the help of Eq.d22b as (n(t) 2 ) = — — ^ " * — A(n(f)), or: 

(n(0 a > = £ (2fc + A)(^+ ft- 1)A) j^^^ e-^)^. (73) 
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